/**********************************************************************
 This file estimates the MC models proposed by Dan Millimet in our 
 joint paper for both Uganda and Kenya.  In Uganda we use only the proximity
 specification.  */
 
 
 
 **MC models**
clear all
set more off
cap prog drop mclogit
cap prog drop mclogit_p
cap prog drop mcre
cap prog drop mcre_p
cap log cl
set seed 235935

use "${data}\UgandaKenya20162019Replicatoin.dta", clear


* Defining key variables and locals

replace year = year-2016
gen temp = dist_bordc

replace maize_prod = maize_prod/1000
gen yrmaize = year*maize_prod


* Starting with Uganda analysis 
if 0 {
preserve 

keep if country == "Uganda"


loc y "defor10"

loc x0 "year yrmaize"
loc x11 "postban yrmaize"

g z1 = num_scenes
g z2 = num_scenes*maize_prod


*create means for CRE models*
loc mx11 " "
loc mz1 " "
loc mz2 " "


foreach x in `x11' {
	bys id: egen m`x'=mean(`x')
	loc mx11 = "`mx1' m`x'"
}
	foreach x in `z1' {
		bys id: egen m`x'=mean(`x')
		loc mz1 = "`mz1' m`x'"
}
	foreach x in `z2' {
		bys id: egen m`x'=mean(`x')
		loc mz2 = "`mz1' m`x'"
}
	

*-----------UGANDA-------------------------------------------*

constraint 3 [a0]_cons = -1
constraint 4 [a0]z1 = 0
constraint 5 [a0]z2 = 0


loc i = 1

logit `y'  `x0' i.postban##c.proximity `mx11' if gte1010 == 1, vce(cl code2) nolog
loc b=_b[_cons]
loc fromp "`y':_cons=`b'"
foreach x in 1.postban 1.postban#c.proximity proximity `mx11' `x0' {
	loc b=_b[`x'] 
	loc fromp "`fromp' `y':`x'=`b'"
}
loc fromp "`fromp' a0:_cons=-5 a0:z1=0 a1:_cons=-1 a1:z1=.04"

mclogit `y' i.postban##c.proximity `x0' `mx11' if gte1010 == 1, vce(cl code2) mc(z1 z2) diff from(`fromp') tech(nr 5 dfp 10 bhhh 10) nolog
loc b=[`y']_cons
loc fromp "`y':_cons=`b'"
foreach x in 1.postban proximity 1.postban#c.proximity `mx11' `x0' {
	loc b=[`y']`x'
	loc fromp "`fromp' `y':`x'=`b'"
}

foreach q in 0 1 {
	loc b1=[a`q']_cons
	loc b2=[a`q']z1
	loc fromp "`fromp' a`q':_cons=`b1' a`q':z1=`b2'"		
}

mclogit `y'  i.postban##c.proximity `x0' `mx11' if gte1010 == 1, vce(cl code2) mc(z1 z2) diff from(`fromp')  tech(nr 5 dfp 10 bhhh 10) nolog
mat b2=e(b),0
loc ll=e(ll)
loc a0=[a0]_cons
loc a1=[a1]_cons
g G0p=[a0]_cons
g G1p=[a1]_cons
foreach x in `z1' `z2'  {
	su `x'
	loc a0=`a0'+r(mean)*[a0]`x'
	loc a1=`a1'+r(mean)*[a1]`x'
	replace G0p=G0p+`x'*[a0]`x'
	replace G1p=G1p+`x'*[a1]`x'
}
loc a0=normal(`a0')
loc a1=normal(`a1')
replace G0p=normal(G0p)
replace G1p=normal(G1p)

cap drop xb
predict xb, xb
eststo col`i': margins,  dydx(1.postban) at(proximity=(.002 .05 .50 1)) noestimcheck predict(pr) post
eststo col`i', add(G0 `a0')
eststo col`i', add(G1 `a1')
eststo col`i', add(logL `ll')
loc i=`i'+1

loc y defor30

logit `y'  `x0' i.postban##c.proximity `mx11' if gte1030 == 1, vce(cl code2) nolog
loc b=_b[_cons]
loc fromp "`y':_cons=`b'"
foreach x in 1.postban 1.postban#c.proximity proximity `mx11' `x0' {
	loc b=_b[`x'] 
	loc fromp "`fromp' `y':`x'=`b'"
}
loc fromp "`fromp' a0:_cons=0 a0:z1=-.05 a1:_cons=2.5 a1:z1=-.14"

mclogit `y' i.postban##c.proximity `mx11' `x0'  if gte30 == 1, vce(cl code2) mc(z1 z2) diff from(`fromp')  tech(bhhh 2 nr 5 dfp 10 bhhh 10 ) nolog

mat b2=e(b),0
loc ll=e(ll)
loc a0=[a0]_cons
loc a1=[a1]_cons
replace G0p=[a0]_cons
replace G1p=[a1]_cons
foreach x in `z1' `z2'  {
	su `x'
	loc a0=`a0'+r(mean)*[a0]`x'
	loc a1=`a1'+r(mean)*[a1]`x'
	replace G0p=G0p+`x'*[a0]`x'
	replace G1p=G1p+`x'*[a1]`x'
}
loc a0=normal(`a0')
loc a1=normal(`a1')
replace G0p=normal(G0p)
replace G1p=normal(G1p)

cap drop xb
predict xb, xb
eststo col`i': margins,  dydx(1.postban) at(proximity=(.002 .05 .50 1)) noestimcheck predict(pr) post
eststo col`i', add(G0 `a0')
eststo col`i', add(G1 `a1')
eststo col`i', add(logL `ll')
loc i=`i'+1

loc y defor50

logit `y'  `x0' i.postban##c.proximity `mx11' if gte1050 == 1, vce(cl code2) nolog
loc b=_b[_cons]
loc fromp "`y':_cons=`b'"
foreach x in 1.postban 1.postban#c.proximity proximity `mx11' `x0' {
	loc b=_b[`x'] 
	loc fromp "`fromp' `y':`x'=`b'"
}
loc fromp "`fromp' a0:_cons=-5 a0:z1=0 a1:_cons=-1 a1:z1=.04"

mclogit `y' i.postban##c.proximity `x0' `mx11' if gte1050 == 1, vce(cl code2) mc(z1 z2) diff from(`fromp') tech(nr 5 dfp 10 bhhh 10) nolog
mat b2=e(b),0
loc ll=e(ll)
loc a0=[a0]_cons
loc a1=[a1]_cons
replace G0p=[a0]_cons
replace G1p=[a1]_cons
foreach x in `z1' `z2'  {
	su `x'
	loc a0=`a0'+r(mean)*[a0]`x'
	loc a1=`a1'+r(mean)*[a1]`x'
	replace G0p=G0p+`x'*[a0]`x'
	replace G1p=G1p+`x'*[a1]`x'
}
loc a0=normal(`a0')
loc a1=normal(`a1')
replace G0p=normal(G0p)
replace G1p=normal(G1p)

cap drop xb
predict xb, xb
eststo col`i': margins,  dydx(1.postban) at(proximity=(.002 .05 .50 1)) noestimcheck predict(pr) post
eststo col`i', add(G0 `a0')
eststo col`i', add(G1 `a1')
eststo col`i', add(logL `ll')
loc i=`i'+1


esttab col*  ///
	using "${tables}\UgandaMC.tex", ///
	se nonotes  style(tex)  b(%12.3f) se(%12.3f)  noobs ///
	starlevels(* 0.10 ** 0.05 *** 0.01) label mlabels("" "" "") ///
	 nonumbers  replace  fragment  ///
	 scalar(G0 G1 alpha logL) prehead({ \begin{tabular}{lcc}    ///
	 \hline \hline ///
	 & \multicolumn{1}{c}{10\% threshold} & \multicolumn{1}{c}{30\% threshold} & \multicolumn{1}{c}{50\% threshold}  \\  ///
	 & (1) & (2) (3)  \\ ) ///
	 prefoot( \\ ) ///
		postfoot( \\ \hline \hline \end{tabular} }  ///
		\begin{tablenotes}[para,flushleft] \footnotesize{Column headers indicate sample used. ///
		Standard errors are clustered at the district level. ///
		Marginal effects evaluated at sample means are displayed. $G_{0}$ and $G_{1}$ are the probability of a false positive and negative, respectively, ///
		evaluated at sample means. Misclassification rates depend on the number of L7 cloud-free scenes and its ///
		interaction with cell level maize productivity. A time trend and trend interacted with maize productivity are included in all models. * p $<$.10, ** p$<$ .05, *** p$<$.01.} \end{tablenotes} )

restore
}

*-----------KENYA-------------------------------------------*

if 1{
	preserve
	keep if country == "Kenya"

constraint 3 [a0]_cons = -1
constraint 4 [a0]z1 = 0
constraint 5 [a0]z2 = 0

loc x0 "year yrmaize"
 
*create means for CRE models*

g z1 = num_scenes
g z2 = num_scenes*maize_prod


*create means for CRE models*
loc mx11 " "
loc mz1 " "
loc mz2 " "


foreach x in yrmaize {
	bys id: egen m`x'=mean(`x')
}



loc y defor30
loc i = 1
keep  if gte1030 == 1

foreach a in 0.10 0.20 {
	loc h=`a'*100
	loc aa = ln(`a')
	constraint 1 /lnalpha = `aa'
	constraint 2 [ln_a2]_cons = `aa'
	constraint 3 [a0]_cons = -3
	constraint 4 [a0]z1 = 0
	constraint 5 [a0]z2 = 0
	
	
	scobit `y' i.postban##i.pa `x0' myrmaize, vce(cl code2) constraint(1) nolog
	loc b=[`y']_cons
	loc froms "`y':_cons=`b'"
	foreach x in `x0' 1.postban 1.pa 1.postban#1.pa  myrmaize {
		loc b=[`y']`x'
		loc froms "`froms' `y':`x'=`b'"
	}	
	loc froms "`froms' a1:_cons=0.65 a1:z1=.02"

	mcre `y' i.postban##i.pa `x0' myrmaize, vce(cl code2) mc(z1 z2) diff from(`froms') constraint(1 2 3 4) tech(nr 5 dfp 30 bhhh 10 ) nolog
	loc b=[`y']_cons
	loc froms "`y':_cons=`b'"
	foreach x in `x0' 1.postban 1.pa 1.postban#1.pa  myrmaize {
		loc b=[`y']`x'
		loc froms "`froms' `y':`x'=`b'"
	}	
	loc b1=[a1]_cons
	loc b2=[a1]z1
	loc froms "`froms' a1:_cons=`b1' a1:z1=`b2'"	
	
}

logit `y'  i.postban##i.pa `x0' myrmaize if gte1030 == 1, vce(cl code2)  nolog
loc b=_b[_cons]
loc fromp "`y':_cons=`b'"
foreach x in `x0' 1.postban 1.pa 1.postban#1.pa  {
	loc b=_b[`x'] 
	loc fromp "`fromp' `y':`x'=`b'"
}
loc fromp "`fromp' a0:_cons=0 a0:z1=0 a1:_cons=-2 a1:z1=.04"

mclogit `y'  i.postban##i.pa  `x0' , vce(cl code2) mc(z1 ) diff from(`fromp') tech(nr 5 dfp 10 bhhh 10) nolog


loc b=[`y']_cons
loc fromp "`y':_cons=`b'"
foreach x in `x0' 1.postban  {
	loc b=[`y']`x'
	loc fromp "`fromp' `y':`x'=`b'"
}

foreach q in 0 1 {
	loc b1=[a`q']_cons
	loc b2=[a`q']z1
	loc fromp "`fromp' a`q':_cons=`b1' a`q':z1=`b2'"		
}
mclogit `y'  i.postban##i.pa `x0' myrmaize if gte1030 == 1, vce(cl code2) mc(z1 z2) diff from(`fromp') tech(nr 5 dfp 10 bhhh 10) nolog



mat b2=e(b),0
loc ll=e(ll)
loc a0=[a0]_cons
loc a1=[a1]_cons
g G0p=[a0]_cons
g G1p=[a1]_cons
foreach x in `z1'  {
	su `x'
	loc a0=`a0'+r(mean)*[a0]`x'
	loc a1=`a1'+r(mean)*[a1]`x'
	replace G0p=G0p+`x'*[a0]`x'
	replace G1p=G1p+`x'*[a1]`x'
}
loc a0=normal(`a0')
loc a1=normal(`a1')
replace G0p=normal(G0p)
replace G1p=normal(G1p)

cap drop xb
predict xb, xb
eststo col`i': margins,  dydx(1.postban) at(pa=(0 1)) noestimcheck predict(pr) post
eststo col`i', add(G0 `a0')
eststo col`i', add(G1 `a1')
eststo col`i', add(logL `ll')
loc i=`i'+1
*/


loc y defor30
keep if gte1030 == 1

logit `y'  i.postban##i.pa `x0' myrmaize if gte1030 == 1, vce(cl code2) nolog
loc b=_b[_cons]
loc fromp "`y':_cons=`b'"
foreach x in 1.postban 1.postban#1.pa myrmaize `x0' {
	loc b=_b[`x'] 
	loc fromp "`fromp' `y':`x'=`b'"
}
loc fromp "`fromp' a0:_cons=0 a0:z1=0 a1:_cons=-1 a1:z1=.04"

mclogit `y'  i.postban##i.pa `x0' myrmaize  if gte1030 == 1, vce(cl code2) mc(z1 z2) diff from(`fromp') tech(nr 5 dfp 10 bhhh 10) nolog
loc b=[`y']_cons
loc fromp "`y':_cons=`b'"
foreach x in 1.postban 1.postban#1.pa myrmaize `x0' {
	loc b=[`y']`x'
	loc fromp "`fromp' `y':`x'=`b'"
}

foreach q in 0 1 {
	loc b1=[a`q']_cons
	loc b2=[a`q']z1
	loc fromp "`fromp' a`q':_cons=-2 a`q':z1=`b2'"		
}


mclogit `y'  i.postban##i.pa `x0' myrmaize if gte1030 == 1, vce(cl code2) mc(z1 z2) diff from(`fromp') tech(nr 5 dfp 10 bhhh 10) nolog
mat b2=e(b),0
loc ll=e(ll)
loc a0=[a0]_cons
loc a1=[a1]_cons
replace G0p=[a0]_cons
replace G1p=[a1]_cons
foreach x in `z1' `z2'  {
	su `x'
	loc a0=`a0'+r(mean)*[a0]`x'
	loc a1=`a1'+r(mean)*[a1]`x'
	replace G0p=G0p+`x'*[a0]`x'
	replace G1p=G1p+`x'*[a1]`x'
}
loc a0=normal(`a0')
loc a1=normal(`a1')
replace G0p=normal(G0p)
replace G1p=normal(G1p)

cap drop xb
predict xb, xb
eststo col`i': margins,  dydx(1.postban) at(pa=(0 1)) noestimcheck predict(pr) post
eststo col`i', add(G0 `a0')
eststo col`i', add(G1 `a1')
eststo col`i', add(logL `ll')
loc i=`i'+1

loc y defor50
keep if gte1050 == 1

logit `y' i.postban##i.pa `x0' myrmaize if gte1050 == 1, vce(cl code2) nolog
loc b=_b[_cons]
loc fromp "`y':_cons=`b'"
foreach x in 1.postban 1.postban#1.pa myrmaize `x0' {
	loc b=_b[`x'] 
	loc fromp "`fromp' `y':`x'=`b'"
}
loc fromp "`fromp' a0:_cons=-1 a0:z1=-.06 a0:z2=0 a1:_cons=2 a1:z1=-.05 a1:z2=0"

mclogit `y' i.postban##i.pa `x0' myrmaize  if gte1050 == 1, vce(cl code2) mc(z1 z2) diff from(`fromp') tech(nr 5 dfp 10 bhhh 10) nolog
loc b=[`y']_cons
loc fromp "`y':_cons=`b'"
foreach x in 1.postban 1.postban#1.pa myrmaize `x0' {
	loc b=[`y']`x'
	loc fromp "`fromp' `y':`x'=`b'"
}

foreach q in 0 1 {
	loc b1=[a`q']_cons
	loc b2=[a`q']z1
	loc fromp "`fromp' a`q':_cons=`b1' a`q':z1=`b2'"		
}


mclogit `y'  i.postban##i.pa `x0' myrmaize  if gte1050 == 1, vce(cl code2) mc(z1 z2) diff from(`fromp') tech(nr 5 dfp 10 bhhh 10) nolog
mat b2=e(b),0
loc ll=e(ll)
loc a0=[a0]_cons
loc a1=[a1]_cons
replace G0p=[a0]_cons
replace G1p=[a1]_cons
foreach x in `z1' `z2'  {
	su `x'
	loc a0=`a0'+r(mean)*[a0]`x'
	loc a1=`a1'+r(mean)*[a1]`x'
	replace G0p=G0p+`x'*[a0]`x'
	replace G1p=G1p+`x'*[a1]`x'
}
loc a0=normal(`a0')
loc a1=normal(`a1')
replace G0p=normal(G0p)
replace G1p=normal(G1p)

cap drop xb
predict xb, xb
eststo col`i': margins,  dydx(1.postban) at(pa=(0  1)) noestimcheck predict(pr) post
eststo col`i', add(G0 `a0')
eststo col`i', add(G1 `a1')
eststo col`i', add(logL `ll')
loc i=`i'+1


esttab col*  ///
	using "${tables}\KenyaMC.tex", ///
	se nonotes  style(tex)  b(%12.3f) se(%12.3f)  noobs ///
	starlevels(* 0.10 ** 0.05 *** 0.01) label mlabels("" "" "" ) ///
	 nonumbers  replace  fragment  ///
	 scalar(G0 G1 alpha logL) prehead({ \begin{tabular}{lccc}    ///
	 \hline \hline ///
	  \multicolumn{1}{c}{10\% threshold} & & \multicolumn{1}{c}{30\% threshold} & \multicolumn{1}{c}{50\% threshold}  \\  ///
	 & (1) & (2) & (3)  \\ ) ///
	 prefoot( \\ ) ///
		postfoot(Full MC covars & no & no & yes \\ \hline \hline \end{tabular} }  ///
		\begin{tablenotes}[para,flushleft] \footnotesize{Column headers indicate estimator. ///
		Standard errors are clustered at the district level. ///
		Marginal effects evaluated at sample means are displayed. $G_{0}$ and $G_{1}$ are the probability of a false positive and negative, respectively, ///
		evaluated at sample means. The misclassification rates depend on the number of L7 cloud-free scenes and its ///
		interaction with cell level maize productivity. A time trend and trend interacted with maize productivity is included in all models. * p $<$.10, ** p$<$ .05, *** p$<$.01.} \end{tablenotes} )
restore

}
	
	
	
